## This is R code for "The Demand for Insurance: Incorporating the Severity of
## Losing Office into the Insurance Model of Judicial Independence"

## 0. Preparation ------------------------------------------------------------
rm(list=ls())

## install.packages("pacman")
pacman::p_load(
  tidyverse, this.path, simcf, haven, runner, arm
)

## ----------------------------------------------------------------------------
## set working directory
setwd(this.path::here(..=1))

## 1. Read the main variables  -------------------------------------------------
dat <- read.csv("Data/dat_2024.csv") 

## 2. Get leader exit from archigos  -------------------------------------------
## To create the punishment variables
archigos <- read_dta("Data/Archigos_4.1_stata14.dta") %>%
  mutate(startyear=year(startdate), endyear=year(enddate)) %>%
  dplyr::select(c(ccode, leader, startyear, endyear, exit, exitcode, posttenurefate)) 
## Turn archigos into country-year_leader format  
archigos_cyl <- archigos %>%  tidyr::expand(nesting(ccode, leader), year=seq(1840,2015)) %>% ## tweak here for start and end year
  left_join(., archigos, by=c("ccode","leader")) %>%
  filter(year >= startyear & year <= endyear) %>% arrange(ccode, year, leader) %>%
  mutate(exitoccur = (year==endyear)&(exit!="Still in Office"), 
         punish_death = (posttenurefate=="Death")&exitoccur, 
         punish_exile = (posttenurefate=="Exile")&exitoccur, 
         punish_imprisonment = (posttenurefate=="Imprisonment")&exitoccur)
## Note, there could be multiple leaders and leader exits in a given year
## So, summarise into country-year format first
archigos_cy <- archigos_cyl %>% 
  filter(year>=1945) %>% # use only post WWII data
  group_by(ccode, year) %>%
  summarise(exit = sum(exitoccur, na.rm = TRUE),
            punish_death = sum(punish_death, na.rm = TRUE),
            punish_exile = sum(punish_exile, na.rm = TRUE),
            punish_imprisonment = sum(punish_imprisonment, na.rm = TRUE)) %>% 
  group_by(ccode) %>%
  mutate(exit_cumsum = cumsum(exit),
         exit_roll10 = runner(exit, sum, k = 10),
         exit_roll20 = runner(exit, sum, k = 20),
         punish_death_cumsum = cumsum(punish_death),
         punish_death_roll10 = runner(punish_death, sum, k = 10),
         punish_death_roll20 = runner(punish_death, sum, k = 20),
         punish_exile_cumsum = cumsum(punish_exile),
         punish_exile_roll10 = runner(punish_exile, sum, k = 10),
         punish_exile_roll20 = runner(punish_exile, sum, k = 20),
         punish_imprisonment_cumsum = cumsum(punish_imprisonment),
         punish_imprisonment_roll10 = runner(punish_imprisonment, sum, k = 10),
         punish_imprisonment_roll20 = runner(punish_imprisonment, sum, k = 20)) %>%
  ungroup() 
## Create percentage of punishment
## weight low: 0.5, 0.25 (imprisonment, exile)
## weight high: 0.75, 0.5
wt_low <- c(0.5, 0.25)
wt_high <- c(0.75, 0.5)
archigos_cy_punish <- archigos_cy %>%
  mutate(punish_percent_cumsum = ifelse(exit_cumsum==0, 0, 
                                        (punish_death_cumsum+punish_exile_cumsum+punish_imprisonment_cumsum)/exit_cumsum),
         punish_percent_roll10 = ifelse(exit_roll10==0, 0, 
                                        (punish_death_roll10+punish_exile_roll10+punish_imprisonment_roll10)/exit_roll10),
         punish_percent_roll20 = ifelse(exit_roll20==0, 0, 
                                        (punish_death_roll20+punish_exile_roll20+punish_imprisonment_roll20)/exit_roll20),
         punish_percent_cumsum_wt_low = ifelse(exit_cumsum==0, 0, 
                                        (punish_death_cumsum+wt_low[2]*punish_exile_cumsum+wt_low[1]*punish_imprisonment_cumsum)/exit_cumsum),
         punish_percent_roll10_wt_low = ifelse(exit_roll10==0, 0, 
                                        (punish_death_roll10+wt_low[2]*punish_exile_roll10+wt_low[1]*punish_imprisonment_roll10)/exit_roll10),
         punish_percent_roll20_wt_low = ifelse(exit_roll20==0, 0, 
                                        (punish_death_roll20+wt_low[2]*punish_exile_roll20+wt_low[1]*punish_imprisonment_roll20)/exit_roll20),
         punish_percent_cumsum_wt_high = ifelse(exit_cumsum==0, 0, 
                                        (punish_death_cumsum+wt_high[2]*punish_exile_cumsum+wt_high[1]*punish_imprisonment_cumsum)/exit_cumsum),
         punish_percent_roll10_wt_high = ifelse(exit_roll10==0, 0, 
                                        (punish_death_roll10+wt_high[2]*punish_exile_roll10+wt_high[1]*punish_imprisonment_roll10)/exit_roll10),
         punish_percent_roll20_wt_high = ifelse(exit_roll20==0, 0, 
                                        (punish_death_roll20+wt_high[2]*punish_exile_roll20+wt_high[1]*punish_imprisonment_roll20)/exit_roll20)
         )

## 3. Merge all the variables together and create logs and lags-----------------
dat <- dat %>% 
  left_join(archigos_cy_punish, by=c("ccode","year")) %>%
  mutate(
    logdemdur = log(e_democracy_duration),
    gdpS = loggdp_lag1 / (2 * sd(loggdp_lag1, na.rm = TRUE)),
    demdurS = logdemdur / (2 * sd(logdemdur, na.rm = TRUE)),
    ross_oil_value_2000_lag1 = lagpanel(ross_oil_value_2000, country_name, year, 1), 
    aid_crsio_lag1 = lagpanel(aid_crsio, country_name, year, 1),
    aid_crsc_lag1 = lagpanel(aid_crsc, country_name, year, 1),   
    e_reserves_billions_lag1 = lagpanel(e_reserves_billions, country_name, year, 1))

## 4. Create neighbor's mean punishment rate -----------------------------------
## read contiguity data 
contdird <- read_csv("Data/contdird.csv") %>%
  dplyr::select(c(state1no, state2no, year))
## Function to calculate neighbor averages
neighbormean <- function(dataset=dat, var="punish_percent_cumsum", 
                         ccodename="ccode", yearname="year",
                         country=2, time=2000){
  neighboraverage = contdird %>% filter(state1no==country, year==time) %>%
    left_join(dataset %>% dplyr::select(one_of(c(ccodename, yearname, var))), 
              by=c("state2no"=ccodename, "year"=yearname)) %>%
    dplyr::select(one_of(var)) %>% unlist() %>% mean(., na.rm=T) # tweak this function to get other summaries
  return(neighboraverage)
}
## generating means for each variables
dat <- dat %>% group_by(year, ccode) %>%
  mutate(neighbormean_punish_percent_cumsum = neighbormean(dataset=dat, var="punish_percent_cumsum", 
                                                     country=ccode, time=year),
         neighbormean_punish_punish_percent_roll10 = neighbormean(dataset=dat, var="punish_percent_roll10", 
                                                           country=ccode, time=year),
         neighbormean_punish_punish_percent_roll20 = neighbormean(dataset=dat, var="punish_percent_roll20", 
                                                                 country=ccode, time=year),
         neighbormean_punish_percent_cumsum_wt_low = neighbormean(dataset=dat, var="punish_percent_cumsum_wt_low", 
                                                           country=ccode, time=year),
         neighbormean_punish_punish_percent_roll10_wt_low = neighbormean(dataset=dat, var="punish_percent_roll10_wt_low", 
                                                                  country=ccode, time=year),
         neighbormean_punish_punish_percent_roll20_wt_low = neighbormean(dataset=dat, var="punish_percent_roll20_wt_low", 
                                                                  country=ccode, time=year),
         neighbormean_punish_percent_cumsum_wt_high = neighbormean(dataset=dat, var="punish_percent_cumsum_wt_high", 
                                                                  country=ccode, time=year),
         neighbormean_punish_punish_percent_roll10_wt_high = neighbormean(dataset=dat, var="punish_percent_roll10_wt_high", 
                                                                         country=ccode, time=year),
         neighbormean_punish_punish_percent_roll20_wt_high = neighbormean(dataset=dat, var="punish_percent_roll20_wt_high", 
                                                                         country=ccode, time=year)) %>%
  ungroup() %>%
  mutate(threat_cumsum = (punish_percent_cumsum + 1) * e_h_polcon3_lag1,
         threatS_cumsum = threat_cumsum / (2 * sd(threat_cumsum, na.rm = TRUE)),
         threat_roll10 = (punish_percent_roll10 + 1) * e_h_polcon3_lag1,
         threatS_roll10 = threat_roll10 / (2 * sd(threat_roll10, na.rm = TRUE)),
         threat_roll20 = (punish_percent_roll20 + 1) * e_h_polcon3_lag1,
         threatS_roll20 = threat_roll20 / (2 * sd(threat_roll20, na.rm = TRUE)),
         threat_cumsum_wt_low  = (punish_percent_cumsum_wt_low  + 1) * e_h_polcon3_lag1,
         threatS_cumsum_wt_low  = threat_cumsum_wt_low / (2 * sd(threat_cumsum_wt_low, na.rm = TRUE)),
         threat_roll10_wt_low = (punish_percent_roll10_wt_low + 1) * e_h_polcon3_lag1,
         threatS_roll10_wt_low = threat_roll10_wt_low / (2 * sd(threat_roll10_wt_low, na.rm = TRUE)),
         threat_roll20_wt_low = (punish_percent_roll20_wt_low + 1) * e_h_polcon3_lag1,
         threatS_roll20_wt_low = threat_roll20_wt_low / (2 * sd(threat_roll20_wt_low, na.rm = TRUE)),
         threat_cumsum_wt_high  = (punish_percent_cumsum_wt_high  + 1) * e_h_polcon3_lag1,
         threatS_cumsum_wt_high  = threat_cumsum_wt_high / (2 * sd(threat_cumsum_wt_high, na.rm = TRUE)),
         threat_roll10_wt_high = (punish_percent_roll10_wt_high + 1) * e_h_polcon3_lag1,
         threatS_roll10_wt_high = threat_roll10_wt_high / (2 * sd(threat_roll10_wt_high, na.rm = TRUE)),
         threat_roll20_wt_high = (punish_percent_roll20_wt_high + 1) * e_h_polcon3_lag1,
         threatS_roll20_wt_high = threat_roll20_wt_high / (2 * sd(threat_roll20_wt_high, na.rm = TRUE)))
  
save(dat, file = "Data/dat.RData")